/*
  File:    DiffusionSolver.cpp
  Purpose: Provides a fully nonlinear diffusional solver 
           for a single time iteration. In particular, it 
		  solves the problem 
		   
		  u_t = div( DiffusionConstant(u,x,y) . grad(u) ) + Sources(u,x,y)
		   
		  with boundary conditions on the left, right, bottom, 
		  and top, of type:
		    
		  'D': Dirichlet
		  'N': Neuman
		  'E': Extrapolation/continuity

           Assumed:
            Can evaluate DiffusionConstant, Source, and 
		   boundary conditions at any point (x,y) where 
		   appropriate. 
    
            DiffusionConstant and Source can be nonlinear 
		   in u. 
   
           Result: 
            The solution u(x,y,t) at t+dt.   
		   
  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cmath>

using namespace std;

#include "Matrix.h"
#include "ImprovedMath.h"
#include "band.h"

extern int Xnodes;
extern int Ynodes;
extern Matrix Phi;

#include "stopwatch.h"

Matrix NonlinearDiffusionIteration( 
   Matrix& X, Matrix& Y, 
   Matrix& PreviousTimeStep , double dt, 
   double (*LeftBC)(double),   char LeftBCtype,  
   double (*RightBC)(double),  char RightBCtype,     
   double (*BottomBC)(double), char BottomBCtype,   
   double (*TopBC)(double),    char TopBCtype,  
   double (*DiffusionConstant)(double,double,double,int,int), 
   double (*Source)(double,double,double,int,int), 
   int MaxIterations, double Tolerance )
{
 Matrix NonzeroBandsList;
 NonzeroBandsList.SetSize(1,9);
	extern Matrix Debuger;
 int Yoffset = PreviousTimeStep.TellRows(); 
 int Size = (PreviousTimeStep.TellRows()) * (PreviousTimeStep.TellCols());
 
 *NonzeroBandsList(0) = 0;
 *NonzeroBandsList(1) = Yoffset;
 *NonzeroBandsList(2) = 2*Yoffset-2;
 *NonzeroBandsList(3) = 2*Yoffset-1;
 *NonzeroBandsList(4) = 2*Yoffset;
 *NonzeroBandsList(5) = 2*Yoffset+1;
 *NonzeroBandsList(6) = 2*Yoffset+2;
 *NonzeroBandsList(7) = 3*Yoffset;
 *NonzeroBandsList(8) = 4*Yoffset;
 
 BandMatrix Coefficients( Size,Size, 2*Yoffset,2*Yoffset, NonzeroBandsList);
 Matrix RHS = zeros(1,Size);
 
 double dx = *X(1) - *X(0);
 double dy = *Y(1) - *Y(0);
 
 double ConstX = dt/square(dx);
 double ConstY = dt/square(dy);
 
 int i,j,N;
 double x,y,u;
 
 // x-direction first;
 
 for( i=0 ; i < PreviousTimeStep.TellRows() ; i++ )
 {
  j=0;
  N = Yoffset*j + i;
  x = *X(i);
  y = *Y(j);
  if( BottomBCtype == 'D' || BottomBCtype == 'd' )
  {
   *Coefficients.eNoUpdate(N,N) += ConstX;
   *RHS(N) += ConstX*BottomBC( x );
  }
  if( BottomBCtype == 'N' || BottomBCtype == 'n' )
  {
   *Coefficients.eNoUpdate(N,N)         -= ConstY;
   *Coefficients.eNoUpdate(N,N+Yoffset) += ConstY;
   *RHS(N) += BottomBC( x ) * ConstY * dy; 
  }
  if( BottomBCtype == 'E' || BottomBCtype == 'e' )
  {
   *Coefficients.eNoUpdate(N,N)           += ConstY;
   *Coefficients.eNoUpdate(N,N+Yoffset)   -= 2.0*ConstY;
   *Coefficients.eNoUpdate(N,N+2*Yoffset) += ConstY;
  }

  j = PreviousTimeStep.TellCols()-1; 
  N = Yoffset*j + i;
  x = *X(i);
  y = *Y(j);
  if( TopBCtype == 'D' )
  {
   *Coefficients.eNoUpdate(N,N) += ConstX;
   *RHS(N) += ConstX*TopBC( x ); 
  }
  if( TopBCtype == 'N' )
  {
   *Coefficients.eNoUpdate(N,N)         += ConstY;
   *Coefficients.eNoUpdate(N,N-Yoffset) -= ConstY;
   *RHS(N) += TopBC( x ) * ConstY * dy; 
  }
  if( TopBCtype == 'E' || TopBCtype == 'e' )
  {
   *Coefficients.eNoUpdate(N,N)           += ConstY;
   *Coefficients.eNoUpdate(N,N-Yoffset)   -= 2.0*ConstY;
   *Coefficients.eNoUpdate(N,N-2*Yoffset) += ConstY;
  }
 }
 
 for( j=1 ; j < PreviousTimeStep.TellCols()-1 ; j++)
 {
  for( i=1 ; i < PreviousTimeStep.TellRows()-1 ; i++ )
  {
   N = Yoffset*j + i;
   x = *X(i);
   y = *Y(j);
   u = *PreviousTimeStep(i,j); 

   double DiffusionLeft  = DiffusionConstant(u,x-0.5*dx,y,i,j);
   double DiffusionRight = DiffusionConstant(u,x+0.5*dx,y,i,j);

   *Coefficients.e(N,N-1) += DiffusionLeft*ConstX;
   *Coefficients.e(N,N) -= (DiffusionLeft+DiffusionRight)*ConstX;
   *Coefficients.e(N,N+1) += DiffusionRight*ConstX;
   
   // finish the backwards time derivative!
   
   *Coefficients.e(N,N) -= 1.0;
   *RHS(N) -= u; 

   // do the source term 
   *RHS(N) -= dt*Source(u,x,y,i,j); 
	  *Debuger(i,j) = Source(u,x,y,i,j); 

	 	
  } 
 }
 
 // y-direction next;
 
// for( j=0 ; j < PreviousTimeStep.TellCols() ; j++ )
 for( j=1 ; j < PreviousTimeStep.TellCols()-1 ; j++ )
 {
  i=0;
  N = Yoffset*j + i;
  x = *X(i);
  y = *Y(j);
  if( LeftBCtype == 'D' )
  {
   *Coefficients.eNoUpdate(N,N) += ConstX;
   *RHS(N) += ConstX*LeftBC( y );
  }
  if( LeftBCtype == 'N' )
  {
   *Coefficients.eNoUpdate(N,N)   -= ConstX;
   *Coefficients.eNoUpdate(N,N+1) += ConstX;
   *RHS(N) += LeftBC( y ) * ConstX* dx; 
  }
  if( LeftBCtype == 'E' )
  {
   *Coefficients.eNoUpdate(N,N)   += ConstX;
   *Coefficients.eNoUpdate(N,N+1) -= 2.0*ConstX;
   *Coefficients.eNoUpdate(N,N+2) += ConstX;
  }

  i = PreviousTimeStep.TellRows()-1; 
  N = Yoffset*j + i;
  x = *X(i);
  y = *Y(j);
  if( RightBCtype == 'D' )
  {
   *Coefficients.eNoUpdate(N,N) += ConstX;
   *RHS(N) += ConstX*RightBC( y ); 
  }
  if( RightBCtype == 'N' )
  {
   *Coefficients.eNoUpdate(N,N)   += ConstX;
   *Coefficients.eNoUpdate(N,N-1) -= ConstX;
   *RHS(N) += RightBC( y ) * ConstX * dx; 
  }
  if( RightBCtype == 'E' )
  {
   *Coefficients.eNoUpdate(N,N)   += ConstX;
   *Coefficients.eNoUpdate(N,N-1) -= 2.0*ConstX;
   *Coefficients.eNoUpdate(N,N-2) += ConstX;
  }
  

 }
  
 for( i=1 ; i < PreviousTimeStep.TellRows()-1 ; i++ )
 {
  for( j=1 ; j < PreviousTimeStep.TellCols()-1 ; j++)
  {
   N = Yoffset*j + i;
   x = *X(i);
   y = *Y(j);
   u = *PreviousTimeStep(i,j); 

   double DiffusionLeft  = DiffusionConstant(u,x,y-0.5*dy,i,j);
   double DiffusionRight = DiffusionConstant(u,x,y+0.5*dy,i,j);
   
   *Coefficients.e(N,N-Yoffset) += DiffusionLeft*ConstY;
   *Coefficients.e(N,N) -= (DiffusionLeft+DiffusionRight)*ConstY;
   *Coefficients.e(N,N+Yoffset) += DiffusionRight*ConstY;
  } 
 }

 // prepare to invert the system!

 Matrix InitialGuess = zeros( 1, RHS.TellCols() );
 for( j=0 ; j < PreviousTimeStep.TellCols() ; j++ )
 {
  for( i=0 ; i < PreviousTimeStep.TellRows() ; i++ )
  {
   N = Yoffset*j + i;
   *InitialGuess(N) = *PreviousTimeStep(i,j);
  }
 }

 // invert the linear system
 
 int OrderL = 5;
/*
 double Tolerance = IntPow( dx , 4 );
 int MaxIterations = 100;
*/ 
 Coefficients.bandlinsolveBiCGStabOfL( 
         &RHS, &InitialGuess, OrderL, Tolerance,MaxIterations);
 
 // update the guess!

 N = 0; 
 Matrix Output = zeros( PreviousTimeStep.TellRows() , PreviousTimeStep.TellCols() );
 
 for( j=0 ; j < PreviousTimeStep.TellCols() ; j++ )
 {
  for( i=0 ; i < PreviousTimeStep.TellRows() ; i++ )
  {
   *Output(i,j) = *InitialGuess(N); N++;
  }
 }
 
 return Output;
}   
Matrix NonlinearConvectionIteration( 
								   Matrix& X, Matrix& Y, 
								   Matrix& PreviousTimeStep , Matrix& AnotherSource ,double dt, 
								   double (*LeftBC)(double),   char LeftBCtype,  
								   double (*RightBC)(double),  char RightBCtype,     
								   double (*BottomBC)(double), char BottomBCtype,   
								   double (*TopBC)(double),    char TopBCtype,  
								   double (*DiffusionConstant)(double,double,double,int,int), 
								   double (*Source)(double,int,int), 
								   int MaxIterations, double Tolerance )
{
	extern Matrix Pressure;
	extern Matrix PressureGradientX;
	extern Matrix PressureGradientY;
	extern double AdvectionConstant;
	extern Matrix BloodVesselIndicator;
	Matrix NonzeroBandsList;
	NonzeroBandsList.SetSize(1,9);
	extern Matrix Debuger;
	int Yoffset = PreviousTimeStep.TellRows(); 
	int Size = (PreviousTimeStep.TellRows()) * (PreviousTimeStep.TellCols());
	
	*NonzeroBandsList(0) = 0;
	*NonzeroBandsList(1) = Yoffset;
	*NonzeroBandsList(2) = 2*Yoffset-2;
	*NonzeroBandsList(3) = 2*Yoffset-1;
	*NonzeroBandsList(4) = 2*Yoffset;
	*NonzeroBandsList(5) = 2*Yoffset+1;
	*NonzeroBandsList(6) = 2*Yoffset+2;
	*NonzeroBandsList(7) = 3*Yoffset;
	*NonzeroBandsList(8) = 4*Yoffset;
	
	BandMatrix Coefficients( Size,Size, 2*Yoffset,2*Yoffset, NonzeroBandsList);
	Matrix RHS = zeros(1,Size);
	
	double dx = *X(1) - *X(0);
	double dy = *Y(1) - *Y(0);
	
	double ConstX = dt/square(dx);
	double ConstY = dt/square(dy);
	
	int i,j,N;
	double x,y,u;
	
	// x-direction first;
	i = 0;
	j = 0;
	N = Yoffset*j + i;
	x = *X(i);
	y = *Y(j);
	
	u = *PreviousTimeStep(i,j); 
	
	double DiffusionLeft  = DiffusionConstant(u,x+0.5*dx,y,i,j);
	double DiffusionRight = DiffusionConstant(u,x+0.5*dx,y,i,j);
	*Coefficients.e(N,N+1) += (DiffusionLeft+DiffusionRight)*ConstX;
	*Coefficients.e(N,N) -= ConstX*(DiffusionLeft+DiffusionRight);
	DiffusionLeft  = DiffusionConstant(u,x,y+0.5*dy,i,j);
	DiffusionRight = DiffusionConstant(u,x,y+0.5*dy,i,j);
	*Coefficients.e(N,N+Yoffset) += (DiffusionLeft+DiffusionRight)*ConstY;
	*Coefficients.e(N,N) -= ConstY*(DiffusionLeft+DiffusionRight);
	*Coefficients.e(N,N) -= 1.0;
	
	*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i+1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
	*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i,j+1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
	*RHS(N) -= u;	
	
	if(*BloodVesselIndicator(i,j)>0)
	{
		*RHS(N) -= (*AnotherSource(i,j)-u);
	}
	*RHS(N) -= dt*Source(u,i,j);
	*Debuger(i,j) = *RHS(N);
	
		
	
	i = 0;
	j = Yoffset-1;
	N = Yoffset*j + i;
	x = *X(i);
	y = *Y(j);
	
	u = *PreviousTimeStep(i,j);
	 DiffusionLeft  = DiffusionConstant(u,x+0.5*dx,y,i,j);
	 DiffusionRight = DiffusionConstant(u,x+0.5*dx,y,i,j);
	*Coefficients.e(N,N+1) += (DiffusionLeft+DiffusionRight)*ConstX;
	*Coefficients.e(N,N) -= ConstX*(DiffusionLeft+DiffusionRight);
	DiffusionLeft  = DiffusionConstant(u,x,y-0.5*dy,i,j);
	DiffusionRight = DiffusionConstant(u,x,y-0.5*dy,i,j);
	*Coefficients.e(N,N-Yoffset) += (DiffusionLeft+DiffusionRight)*ConstY;
	*Coefficients.e(N,N) -= ConstY*(DiffusionLeft+DiffusionRight);
	*Coefficients.e(N,N) -= 1.0;
	
	*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i+1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
	*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j-1));
	*RHS(N) -= u;	
	
	if(*BloodVesselIndicator(i,j)>0)
	{
		*RHS(N) -= (*AnotherSource(i,j)-u);
	}
	*RHS(N) -= dt*Source(u,i,j);
	*Debuger(i,j) = *RHS(N);
	
	
	
	i = Yoffset-1;
	j = 0;
	N = Yoffset*j + i;
	x = *X(i);
	y = *Y(j);
	u = *PreviousTimeStep(i,j);
	DiffusionLeft  = DiffusionConstant(u,x-0.5*dx,y,i,j);
    DiffusionRight = DiffusionConstant(u,x-0.5*dx,y,i,j);
	*Coefficients.e(N,N-1) += (DiffusionLeft+DiffusionRight)*ConstX;
	*Coefficients.e(N,N) -= ConstX*(DiffusionLeft+DiffusionRight);
	DiffusionLeft  = DiffusionConstant(u,x,y+0.5*dy,i,j);
	DiffusionRight = DiffusionConstant(u,x,y+0.5*dy,i,j);
	*Coefficients.e(N,N+Yoffset) += (DiffusionLeft+DiffusionRight)*ConstY;
	*Coefficients.e(N,N) -= ConstY*(DiffusionLeft+DiffusionRight);
	*Coefficients.e(N,N) -= 1.0;
	
	*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i-1,j));
	*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i,j+1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
	*RHS(N) -= u;	
	
	if(*BloodVesselIndicator(i,j)>0)
	{
		*RHS(N) -= (*AnotherSource(i,j)-u);
	}
	*RHS(N) -= dt*Source(u,i,j);
	*Debuger(i,j) = *RHS(N);
	
	
	i = Yoffset-1;
	j = Yoffset-1;
	N = Yoffset*j + i;
	x = *X(i);
	y = *Y(j);
	u = *PreviousTimeStep(i,j);
	 DiffusionLeft  = DiffusionConstant(u,x-0.5*dx,y,i,j);
	 DiffusionRight = DiffusionConstant(u,x-0.5*dx,y,i,j);
	*Coefficients.e(N,N-1) += (DiffusionLeft+DiffusionRight)*ConstX;
	*Coefficients.e(N,N) -= ConstX*(DiffusionLeft+DiffusionRight);
	DiffusionLeft  = DiffusionConstant(u,x,y-0.5*dy,i,j);
	DiffusionRight = DiffusionConstant(u,x,y-0.5*dy,i,j);
	*Coefficients.e(N,N-Yoffset) += (DiffusionLeft+DiffusionRight)*ConstY;
	*Coefficients.e(N,N) -= ConstY*(DiffusionLeft+DiffusionRight);
	*Coefficients.e(N,N) -= 1.0;
	
	*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i-1,j));
	*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j-1));
	*RHS(N) -= u;	
	
	if(*BloodVesselIndicator(i,j)>0)
	{
		*RHS(N) -= (*AnotherSource(i,j)-u);
	}
	*RHS(N) -= dt*Source(u,i,j);
	*Debuger(i,j) = *RHS(N);
	
	
	for( i=1 ; i < PreviousTimeStep.TellRows()-1 ; i++ )
	{
		j=0;
		N = Yoffset*j + i;
		x = *X(i);
		y = *Y(j);
		if( BottomBCtype == 'D' || BottomBCtype == 'd' )
		{
			*Coefficients.eNoUpdate(N,N) += ConstX;
			*RHS(N) += ConstX*BottomBC( x );
		}
		if( BottomBCtype == 'N' || BottomBCtype == 'n' )
		{
			//*Coefficients.eNoUpdate(N,N)         -= ConstY;
			//*Coefficients.eNoUpdate(N,N+Yoffset) += ConstY;
			//*RHS(N) += BottomBC( x ) * ConstY * dy; 
			
			
			
			N = Yoffset*j + i;
			x = *X(i);
			y = *Y(j);
			
			u = *PreviousTimeStep(i,j); 
			
			double DiffusionLeft  = DiffusionConstant(u,x-0.5*dx,y,i,j);
			double DiffusionRight = DiffusionConstant(u,x+0.5*dx,y,i,j);
			
			
			//*Coefficients.e(N,N-1) += (DiffusionLeft+0.5*AdvectionConstant*(*Pressure(i-1,j)-*Pressure(i,j)))*ConstX;
			//*Coefficients.e(N,N) -= (DiffusionLeft+DiffusionRight-0.5*AdvectionConstant*(*Pressure(i-1,j)-*Pressure(i,j))-0.5*AdvectionConstant*(*Pressure(i+1,j)-*Pressure(i,j)))*ConstX;
			//*Coefficients.e(N,N+1) += (DiffusionRight+0.5*AdvectionConstant*(*Pressure(i+1,j)-*Pressure(i,j)))*ConstX;
			//*Coefficients.e(N,N-1) += (DiffusionLeft+0.5*dx*AdvectionConstant* *PressureGradientX(i,j))*ConstX;
			//*Coefficients.e(N,N+1) += (DiffusionRight-0.5*dx*AdvectionConstant* *PressureGradientX(i,j))*ConstX;
			*Coefficients.e(N,N-1) += (DiffusionLeft)*ConstX;
			*Coefficients.e(N,N+1) += (DiffusionRight)*ConstX;
			*Coefficients.e(N,N) -= 1.0;
			//*Coefficients.e(N,N) += AdvectionConstant*dt*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
			*Coefficients.e(N,N) -= ConstX*(DiffusionLeft+DiffusionRight);
			
			*RHS(N) -= dt*Source(u,i,j);
			if(*BloodVesselIndicator(i,j)>0)
			{
				*RHS(N) -= (*AnotherSource(i,j)-u);
			}
			
			*RHS(N) -= u;			
			//*RHS(N) -= u*dt*AdvectionConstant*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
			
			//*RHS(N) -= +0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i+1,j);
			//*RHS(N) -= -0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i-1,j);
			
			*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
			
			*RHS(N) += +0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i+1,j);
			*RHS(N) += -0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i-1,j);
			
			
			
			
			
			 DiffusionLeft  = DiffusionConstant(u,x,y+0.5*dy,i,j);
			 DiffusionRight = DiffusionConstant(u,x,y+0.5*dy,i,j);
			//*RHS(N) -= +0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j+1);
			//*RHS(N) -= -0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j-1);
			//*RHS(N) -= u*dt*AdvectionConstant*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
			
			//*RHS(N) += +0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j+1);
			//*RHS(N) += -0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j-1);
			*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i,j+1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
			
			//*Coefficients.e(N,N-Yoffset) += (DiffusionLeft+0.5*AdvectionConstant*(*Pressure(i,j-1)-*Pressure(i,j)))*ConstY;
			//*Coefficients.e(N,N) -= (DiffusionLeft+DiffusionRight-0.5*AdvectionConstant*(*Pressure(i,j-1)-*Pressure(i,j))-0.5*AdvectionConstant*(*Pressure(i,j+1)-*Pressure(i,j)))*ConstY;
			//*Coefficients.e(N,N+Yoffset) += (DiffusionRight+0.5*AdvectionConstant*(*Pressure(i,j+1)-*Pressure(i,j)))*ConstY;
			
			//*Coefficients.e(N,N-Yoffset) += (DiffusionLeft+0.5*dx*AdvectionConstant* *PressureGradientY(i,j))*ConstY;
			//*Coefficients.e(N,N+Yoffset) += (DiffusionRight-0.5*dx*AdvectionConstant* *PressureGradientY(i,j))*ConstY;
			//*Coefficients.e(N,N) += AdvectionConstant*ConstY*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
			//*Coefficients.e(N,N) -= ConstY*(DiffusionLeft+DiffusionRight);
			//*Coefficients.e(N,N-Yoffset) += (DiffusionLeft)*ConstY;
			*Coefficients.e(N,N+Yoffset) += ConstY*(DiffusionLeft+DiffusionRight);
			*Coefficients.e(N,N) -= ConstY*(DiffusionLeft+DiffusionRight);
			*Debuger(i,j) = *RHS(N);
			
			
			
			
			
		}
		if( BottomBCtype == 'E' || BottomBCtype == 'e' )
		{
			*Coefficients.eNoUpdate(N,N)           += ConstY;
			*Coefficients.eNoUpdate(N,N+Yoffset)   -= 2.0*ConstY;
			*Coefficients.eNoUpdate(N,N+2*Yoffset) += ConstY;
		}
		
		
		
		
		j = PreviousTimeStep.TellCols()-1; 
		N = Yoffset*j + i;
		x = *X(i);
		y = *Y(j);
		if( TopBCtype == 'D' )
		{
			*Coefficients.eNoUpdate(N,N) += ConstX;
			*RHS(N) += ConstX*TopBC( x ); 
		}
		if( TopBCtype == 'N' )
		{   //add here
			N = Yoffset*j + i;
			x = *X(i);
			y = *Y(j);
			
			u = *PreviousTimeStep(i,j); 
			
			double DiffusionLeft  = DiffusionConstant(u,x-0.5*dx,y,i,j);
			double DiffusionRight = DiffusionConstant(u,x+0.5*dx,y,i,j);
			
			
			//*Coefficients.e(N,N-1) += (DiffusionLeft+0.5*AdvectionConstant*(*Pressure(i-1,j)-*Pressure(i,j)))*ConstX;
			//*Coefficients.e(N,N) -= (DiffusionLeft+DiffusionRight-0.5*AdvectionConstant*(*Pressure(i-1,j)-*Pressure(i,j))-0.5*AdvectionConstant*(*Pressure(i+1,j)-*Pressure(i,j)))*ConstX;
			//*Coefficients.e(N,N+1) += (DiffusionRight+0.5*AdvectionConstant*(*Pressure(i+1,j)-*Pressure(i,j)))*ConstX;
			//*Coefficients.e(N,N-1) += (DiffusionLeft+0.5*dx*AdvectionConstant* *PressureGradientX(i,j))*ConstX;
			//*Coefficients.e(N,N+1) += (DiffusionRight-0.5*dx*AdvectionConstant* *PressureGradientX(i,j))*ConstX;
			*Coefficients.e(N,N-1) += (DiffusionLeft)*ConstX;
			*Coefficients.e(N,N+1) += (DiffusionRight)*ConstX;
			*Coefficients.e(N,N) -= 1.0;
			//*Coefficients.e(N,N) += AdvectionConstant*dt*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
			*Coefficients.e(N,N) -= ConstX*(DiffusionLeft+DiffusionRight);
			
			*RHS(N) -= dt*Source(u,i,j);
			if(*BloodVesselIndicator(i,j)>0)
			{
				*RHS(N) -= (*AnotherSource(i,j)-u);
			}
			
			*RHS(N) -= u;			
			//*RHS(N) -= u*dt*AdvectionConstant*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
			
			//*RHS(N) -= +0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i+1,j);
			//*RHS(N) -= -0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i-1,j);
			
			*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
			
			*RHS(N) += +0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i+1,j);
			*RHS(N) += -0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i-1,j);
						
			
			
			
			DiffusionLeft  = DiffusionConstant(u,x,y-0.5*dy,i,j);
			DiffusionRight = DiffusionConstant(u,x,y-0.5*dy,i,j);
			//*RHS(N) -= +0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j+1);
			//*RHS(N) -= -0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j-1);
			//*RHS(N) -= u*dt*AdvectionConstant*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
			
			//*RHS(N) += +0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j+1);
			//*RHS(N) += -0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j-1);
			*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j-1));
			
			//*Coefficients.e(N,N-Yoffset) += (DiffusionLeft+0.5*AdvectionConstant*(*Pressure(i,j-1)-*Pressure(i,j)))*ConstY;
			//*Coefficients.e(N,N) -= (DiffusionLeft+DiffusionRight-0.5*AdvectionConstant*(*Pressure(i,j-1)-*Pressure(i,j))-0.5*AdvectionConstant*(*Pressure(i,j+1)-*Pressure(i,j)))*ConstY;
			//*Coefficients.e(N,N+Yoffset) += (DiffusionRight+0.5*AdvectionConstant*(*Pressure(i,j+1)-*Pressure(i,j)))*ConstY;
			
			//*Coefficients.e(N,N-Yoffset) += (DiffusionLeft+0.5*dx*AdvectionConstant* *PressureGradientY(i,j))*ConstY;
			//*Coefficients.e(N,N+Yoffset) += (DiffusionRight-0.5*dx*AdvectionConstant* *PressureGradientY(i,j))*ConstY;
			//*Coefficients.e(N,N) += AdvectionConstant*ConstY*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
			//*Coefficients.e(N,N) -= ConstY*(DiffusionLeft+DiffusionRight);
			*Coefficients.e(N,N-Yoffset) += ConstY*(DiffusionLeft+DiffusionRight);
			//*Coefficients.e(N,N+Yoffset) += ConstY*(DiffusionLeft+DiffusionRight);
			*Coefficients.e(N,N) -= ConstY*(DiffusionLeft+DiffusionRight);
			*Debuger(i,j) = *RHS(N);
			
		}
		if( TopBCtype == 'E' || TopBCtype == 'e' )
		{
			*Coefficients.eNoUpdate(N,N)           += ConstY;
			*Coefficients.eNoUpdate(N,N-Yoffset)   -= 2.0*ConstY;
			*Coefficients.eNoUpdate(N,N-2*Yoffset) += ConstY;
		}
	}
	
	
	for( j=1 ; j < PreviousTimeStep.TellCols()-1 ; j++ )
	{
		i=0;
		N = Yoffset*j + i;
		x = *X(i);
		y = *Y(j);
		if( LeftBCtype == 'D' )
		{
			*Coefficients.eNoUpdate(N,N) += ConstX;
			*RHS(N) += ConstX*LeftBC( y );
		}
		if( LeftBCtype == 'N' )
		{
			//*Coefficients.eNoUpdate(N,N)   -= ConstX;
			//*Coefficients.eNoUpdate(N,N+1) += ConstX;
			//*RHS(N) += LeftBC( y ) * ConstX* dx; 
			
			
			
			N = Yoffset*j + i;
			x = *X(i);
			y = *Y(j);
			
			u = *PreviousTimeStep(i,j); 
			
			double DiffusionLeft  = DiffusionConstant(u,x+0.5*dx,y,i,j);
			double DiffusionRight = DiffusionConstant(u,x+0.5*dx,y,i,j);
			
			
			//*Coefficients.e(N,N-1) += (DiffusionLeft+0.5*AdvectionConstant*(*Pressure(i-1,j)-*Pressure(i,j)))*ConstX;
			//*Coefficients.e(N,N) -= (DiffusionLeft+DiffusionRight-0.5*AdvectionConstant*(*Pressure(i-1,j)-*Pressure(i,j))-0.5*AdvectionConstant*(*Pressure(i+1,j)-*Pressure(i,j)))*ConstX;
			//*Coefficients.e(N,N+1) += (DiffusionRight+0.5*AdvectionConstant*(*Pressure(i+1,j)-*Pressure(i,j)))*ConstX;
			//*Coefficients.e(N,N-1) += (DiffusionLeft+0.5*dx*AdvectionConstant* *PressureGradientX(i,j))*ConstX;
			//*Coefficients.e(N,N+1) += (DiffusionRight-0.5*dx*AdvectionConstant* *PressureGradientX(i,j))*ConstX;
			//*Coefficients.e(N,N-1) += (DiffusionLeft)*ConstX;
			*Coefficients.e(N,N+1) += (DiffusionLeft+DiffusionRight)*ConstX;
			*Coefficients.e(N,N) -= 1.0;
			//*Coefficients.e(N,N) += AdvectionConstant*dt*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
			*Coefficients.e(N,N) -= ConstX*(DiffusionLeft+DiffusionRight);
			
			*RHS(N) -= dt*Source(u,i,j);
			if(*BloodVesselIndicator(i,j)>0)
			{
				*RHS(N) -= (*AnotherSource(i,j)-u);
			}
			
			*RHS(N) -= u;			
			//*RHS(N) -= u*dt*AdvectionConstant*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
			
			//*RHS(N) -= +0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i+1,j);
			//*RHS(N) -= -0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i-1,j);
			
			*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i+1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
			
			//*RHS(N) += +0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i+1,j);
			//*RHS(N) += -0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i-1,j);
			
			
			
			
			
			DiffusionLeft  = DiffusionConstant(u,x,y-0.5*dy,i,j);
			DiffusionRight = DiffusionConstant(u,x,y+0.5*dy,i,j);
			//*RHS(N) -= +0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j+1);
			//*RHS(N) -= -0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j-1);
			//*RHS(N) -= u*dt*AdvectionConstant*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
			
			//*RHS(N) += +0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j+1);
			//*RHS(N) += -0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j-1);
			*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
			
			//*Coefficients.e(N,N-Yoffset) += (DiffusionLeft+0.5*AdvectionConstant*(*Pressure(i,j-1)-*Pressure(i,j)))*ConstY;
			//*Coefficients.e(N,N) -= (DiffusionLeft+DiffusionRight-0.5*AdvectionConstant*(*Pressure(i,j-1)-*Pressure(i,j))-0.5*AdvectionConstant*(*Pressure(i,j+1)-*Pressure(i,j)))*ConstY;
			//*Coefficients.e(N,N+Yoffset) += (DiffusionRight+0.5*AdvectionConstant*(*Pressure(i,j+1)-*Pressure(i,j)))*ConstY;
			
			//*Coefficients.e(N,N-Yoffset) += (DiffusionLeft+0.5*dx*AdvectionConstant* *PressureGradientY(i,j))*ConstY;
			//*Coefficients.e(N,N+Yoffset) += (DiffusionRight-0.5*dx*AdvectionConstant* *PressureGradientY(i,j))*ConstY;
			//*Coefficients.e(N,N) += AdvectionConstant*ConstY*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
			//*Coefficients.e(N,N) -= ConstY*(DiffusionLeft+DiffusionRight);
			*Coefficients.e(N,N-Yoffset) += (DiffusionLeft)*ConstY;
			*Coefficients.e(N,N+Yoffset) += ConstY*(DiffusionRight);
			*Coefficients.e(N,N) -= ConstY*(DiffusionLeft+DiffusionRight);
			*Debuger(i,j) = *RHS(N);
			
		}
		if( LeftBCtype == 'E' )
		{
			*Coefficients.eNoUpdate(N,N)   += ConstX;
			*Coefficients.eNoUpdate(N,N+1) -= 2.0*ConstX;
			*Coefficients.eNoUpdate(N,N+2) += ConstX;
		}
		
		i = PreviousTimeStep.TellRows()-1; 
		N = Yoffset*j + i;
		x = *X(i);
		y = *Y(j);
		if( RightBCtype == 'D' )
		{
			*Coefficients.eNoUpdate(N,N) += ConstX;
			*RHS(N) += ConstX*RightBC( y ); 
		}
		if( RightBCtype == 'N' )
		{
			
			
			
			N = Yoffset*j + i;
			x = *X(i);
			y = *Y(j);
			
			u = *PreviousTimeStep(i,j); 
			
			double DiffusionLeft  = DiffusionConstant(u,x-0.5*dx,y,i,j);
			double DiffusionRight = DiffusionConstant(u,x-0.5*dx,y,i,j);
			
			
			//*Coefficients.e(N,N-1) += (DiffusionLeft+0.5*AdvectionConstant*(*Pressure(i-1,j)-*Pressure(i,j)))*ConstX;
			//*Coefficients.e(N,N) -= (DiffusionLeft+DiffusionRight-0.5*AdvectionConstant*(*Pressure(i-1,j)-*Pressure(i,j))-0.5*AdvectionConstant*(*Pressure(i+1,j)-*Pressure(i,j)))*ConstX;
			//*Coefficients.e(N,N+1) += (DiffusionRight+0.5*AdvectionConstant*(*Pressure(i+1,j)-*Pressure(i,j)))*ConstX;
			//*Coefficients.e(N,N-1) += (DiffusionLeft+0.5*dx*AdvectionConstant* *PressureGradientX(i,j))*ConstX;
			//*Coefficients.e(N,N+1) += (DiffusionRight-0.5*dx*AdvectionConstant* *PressureGradientX(i,j))*ConstX;
			//*Coefficients.e(N,N-1) += (DiffusionLeft)*ConstX;
			*Coefficients.e(N,N-1) += (DiffusionLeft+DiffusionRight)*ConstX;
			*Coefficients.e(N,N) -= 1.0;
			//*Coefficients.e(N,N) += AdvectionConstant*dt*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
			*Coefficients.e(N,N) -= ConstX*(DiffusionLeft+DiffusionRight);
			
			*RHS(N) -= dt*Source(u,i,j);
			if(*BloodVesselIndicator(i,j)>0)
			{
				*RHS(N) -= (*AnotherSource(i,j)-u);
			}
			
			*RHS(N) -= u;			
			//*RHS(N) -= u*dt*AdvectionConstant*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
			
			//*RHS(N) -= +0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i+1,j);
			//*RHS(N) -= -0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i-1,j);
			
			*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i-1,j));
			
			//*RHS(N) += +0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i+1,j);
			//*RHS(N) += -0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i-1,j);
			
			
			
			
			
			DiffusionLeft  = DiffusionConstant(u,x,y-0.5*dy,i,j);
			DiffusionRight = DiffusionConstant(u,x,y+0.5*dy,i,j);
			//*RHS(N) -= +0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j+1);
			//*RHS(N) -= -0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j-1);
			//*RHS(N) -= u*dt*AdvectionConstant*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
			
			//*RHS(N) += +0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j+1);
			//*RHS(N) += -0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j-1);
			*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
			
			//*Coefficients.e(N,N-Yoffset) += (DiffusionLeft+0.5*AdvectionConstant*(*Pressure(i,j-1)-*Pressure(i,j)))*ConstY;
			//*Coefficients.e(N,N) -= (DiffusionLeft+DiffusionRight-0.5*AdvectionConstant*(*Pressure(i,j-1)-*Pressure(i,j))-0.5*AdvectionConstant*(*Pressure(i,j+1)-*Pressure(i,j)))*ConstY;
			//*Coefficients.e(N,N+Yoffset) += (DiffusionRight+0.5*AdvectionConstant*(*Pressure(i,j+1)-*Pressure(i,j)))*ConstY;
			
			//*Coefficients.e(N,N-Yoffset) += (DiffusionLeft+0.5*dx*AdvectionConstant* *PressureGradientY(i,j))*ConstY;
			//*Coefficients.e(N,N+Yoffset) += (DiffusionRight-0.5*dx*AdvectionConstant* *PressureGradientY(i,j))*ConstY;
			//*Coefficients.e(N,N) += AdvectionConstant*ConstY*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
			//*Coefficients.e(N,N) -= ConstY*(DiffusionLeft+DiffusionRight);
			*Coefficients.e(N,N-Yoffset) += (DiffusionLeft)*ConstY;
			*Coefficients.e(N,N+Yoffset) += ConstY*(DiffusionRight);
			*Coefficients.e(N,N) -= ConstY*(DiffusionLeft+DiffusionRight);
			*Debuger(i,j) = *RHS(N);
			
		}
		if( RightBCtype == 'E' )
		{
			*Coefficients.eNoUpdate(N,N)   += ConstX;
			*Coefficients.eNoUpdate(N,N-1) -= 2.0*ConstX;
			*Coefficients.eNoUpdate(N,N-2) += ConstX;
		}
		
		
	}
	
	
	
	for( j=1 ; j < PreviousTimeStep.TellCols()-1 ; j++)
	{
		for( i=1 ; i < PreviousTimeStep.TellRows()-1 ; i++ )
		{
			N = Yoffset*j + i;
			x = *X(i);
			y = *Y(j);
			u = *PreviousTimeStep(i,j); 
			
			double DiffusionLeft  = DiffusionConstant(u,x-0.5*dx,y,i,j);
			double DiffusionRight = DiffusionConstant(u,x+0.5*dx,y,i,j);

			*Coefficients.e(N,N-1) += (DiffusionLeft)*ConstX;
			*Coefficients.e(N,N+1) += (DiffusionRight)*ConstX;
			*Coefficients.e(N,N) -= 1.0;
			*Coefficients.e(N,N) -= ConstX*(DiffusionLeft+DiffusionRight);
			
			*RHS(N) -= dt*Source(u,i,j);
			if(*BloodVesselIndicator(i,j)>0)
			{
			*RHS(N) -= (*AnotherSource(i,j)-u);
			}
			
			*RHS(N) -= u;			
			*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i-1,j)-2* *Pressure(i,j)+ *Pressure(i+1,j));
			
			*RHS(N) += +0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i+1,j);
			*RHS(N) += -0.5*dx*AdvectionConstant* *PressureGradientX(i,j)*ConstX* *PreviousTimeStep(i-1,j);
			 
			
			
			 DiffusionLeft  = DiffusionConstant(u,x,y-0.5*dy,i,j);
			 DiffusionRight = DiffusionConstant(u,x,y+0.5*dy,i,j);

			
			*RHS(N) += +0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j+1);
			*RHS(N) += -0.5*dx*AdvectionConstant* *PressureGradientY(i,j)*ConstY* *PreviousTimeStep(i,j-1);
			*RHS(N) += u*dt*AdvectionConstant*(*Pressure(i,j-1)-2* *Pressure(i,j)+ *Pressure(i,j+1));
			
	
			*Coefficients.e(N,N-Yoffset) += (DiffusionLeft)*ConstY;
			*Coefficients.e(N,N+Yoffset) += (DiffusionRight)*ConstY;
			*Coefficients.e(N,N) -= ConstY*(DiffusionLeft+DiffusionRight);
			*Debuger(i,j) = *RHS(N);
			
		} 
	}
	
	
	// prepare to invert the system!
	
	Matrix InitialGuess = zeros( 1, RHS.TellCols() );
	for( j=0 ; j < PreviousTimeStep.TellCols() ; j++ )
	{
		for( i=0 ; i < PreviousTimeStep.TellRows() ; i++ )
		{
			N = Yoffset*j + i;
			*InitialGuess(N) = *PreviousTimeStep(i,j);
		}
	}
	
	// invert the linear system
	
	int OrderL = 5;
	/*
	 double Tolerance = IntPow( dx , 4 );
	 int MaxIterations = 100;
	 */ 
	Coefficients.bandlinsolveBiCGStabOfL( 
										 &RHS, &InitialGuess, OrderL, Tolerance,MaxIterations);
	
	// update the guess!
	
	N = 0; 
	Matrix Output = zeros( PreviousTimeStep.TellRows() , PreviousTimeStep.TellCols() );
	
	for( j=0 ; j < PreviousTimeStep.TellCols() ; j++ )
	{
		for( i=0 ; i < PreviousTimeStep.TellRows() ; i++ )
		{
			*Output(i,j) = *InitialGuess(N); N++;
		}
	}
	
	return Output;
}   

